Nonequilibrium Invariant Measure under Heat Flow 
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We provide an explicit representation of the nonequilibrium invariant measure for a chain of har- 
monic oscillators with conservative noise in the presence of stationary heat flow. By first determining 
the covariance matrix, we are able to express the measure as the product of Gaussian distributions 
aligned along some collective modes that are spatially localized with power-law tails. Numerical 
studies show that such a representation applies also to a purely deterministic model, the quartic 
Fermi-Pasta-Ulam chain. 
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Characterizing the invariant measure of systems 
steadily kept out of equilibrium is one of the main chal- 
lenges towards the construction of a nonequilibrium sta- 
tistical thermodynamics. Some insight has been gained 
over the years mostly thanks to the solution of stochas- 
tic models In particular, heat conductivity in one- 
dimensional systems has played a prominent role, be- 
cause its understanding would also imply giving a mi- 
croscopic basis to the Fourier's law 0. The first step 
dates back to 1967, when the problem of a harmonic 
chain in contact with two thermal reservoirs at different 
temperatures was solved exactly 0). Unfortunately, the 
integrability of the underlying dynamics leads to several 
unphysical features (e.g. a vanishing temperature gradi- 
ent) |2j. Later, purely stochastic models, where energy 
is assumed to diffuse between neighbouring boxes (the 
oscillators) , have been considered [3, [j| . More recently, 
systems of harmonic oscillators exchanging energy with 
"conservative" noise have been proven to admit a unique 
stationary state with a constant heat flux and a linear 
temperature profile @. Admittedly, the leap from such 
class of models to even the simplest deterministic, non- 
linear ones is still a challenge for the theory [7j ■ 

In this Letter we approach the problem of describing 
the invariant measure in terms of the spectral properties 
of the covariance matrix, i.e. at the level of two-point 
correlators of the relevant dynamical variables. This pro- 
cedure is often referred to in data analysis as principal 
component analysis [§]. It amounts to expressing the 
initial distribution as the product of univariate Gaus- 
sians aligned along the eigenvectors of the covariance 
matrix and whose widths are given by the correspond- 
ing eigenvalues. The decomposition is exact for multi- 
variate Gaussians (linear processes). In order to facil- 
itate the comparison with the Gibbs equilibrium mea- 
sure, it is convenient to express the covariance matrix in 
those variables that make it perfectly diagonal at equi- 
librium. At variance with purely diffusive models, the 
appearance of non-diagonal terms is crucially related to 



the onset of a non-zero heat flux and in particular to the 
existence of anomalous transport properties. Due to the 
almost-diagonal structure, the eigenvectors are localized 
in space. 

To illustrate the approach, we first introduce a stochas- 
tic model, closely related to the one presented in Ref. [j| , 
in which coupled harmonic oscillators interact also by 
random pairwise collisions, each conserving both energy 
and momentum. Due to the linearity of the associated 
master equation, the equations for the two-point correla- 
tors are closed and can be computed exactly without any 
factorization assumption. Another useful property of the 
model is that it encompasses, as a limit case, an anoma- 
lous regime, in which the relevant transport coefficient, 
the thermal conductivity, diverges by virtue of total mo- 
mentum conservation [1 01 ] - This allows comparing the 
cases of normal and anomalous conductivity. In order 
to test the generality of these results, we finally perform 
the same type of analysis by directly simulating a non- 
linear deterministic model. In spite of the unavoidable 
statistical fluctuations, we find encouraging qualitative 
similarities. 

The first model we deal with is a harmonic chain of N 
unit-mass particles, whose displacements and momenta 
are denoted by qi and pi = q\. The chain is in contact 
with stochastic Langevin heat baths at its extrema, 



pi = —kqi + u 2 (qi+i - 2qt + <7i_i) + 
<5i,i(£+ - Api) + Si,N(t,- - Apat), 



(1) 



where £±'s are independent Wiener processes with zero 
mean and variance 2\kgT±, respectively. Fixed bound- 
ary conditions go = 0, (Jn+i = are assumed. Besides 
satisfying dynamics |T]), neighbouring particles are as- 
sumed to randomly collide, and thereby to exchange 
their momenta, with a rate 7. As a result, the phase- 
space probability density P(x,t), that we write in terms 
of the 2iV-dimensional vector whose components 
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(gi, . . . , qw,pi, ■ ■ ■ ,Pn), satisfies the master equation 

(L + L co/ )P . (2) 



dP 
~dt 



The first term accounts for both the deterministic force 
and the coupling with the heat bath and reads 



r) d B 2 P 

a —(x P) + M " 



(3) 



where a M „ and are elements of the 27V x 2 N matrices 
a and d that we write in terms of TV x TV blocks 



a " \to 2 g + kI Ar) ' d 





2\k B T(v + r,s] 



(4) 

with I, being the identity and null matrices respec- 
tively, Tij = Sij(Sn + 5 iN ), Sij = 8 ij {8 i i - 8 iN ) and g i} = 
2§ij - 5 i+ ij - 5 itj+ i. Moreover, we let T = (T+ + T_)/2 
and 7/ = (T + — T_)/T. The collisional term writes 

^coiP = 1^2[P(- ■ -Pi+uPi ■ ■ -)~ p (- ■ -PhPi+i ■ ■ ■)] ( 5 ) 

i 

For k — 7 = the model reduces to the one of Ref. 
where it was showed that the nonequilibrium measure is 
a multivariate Gaussian and all the second moments were 
exactly determined. 

By denoting with (.) the average over P, the covariance 
matrix c^v = (xu,x v ) can be written as four NxN blocks, 



c = 



(6) 



with Uij = (qiqj), % = (piPj), = (qiPj); the symbol f 
denotes the transpose and we also assume = 0. The 
evolution equation is obtained by multiplying both sides 
of Eq. ([2]) by XkXi and thereby integrating, 



c = d ac ca T 



(7) 



The first three terms on the r.h.s. are associated with 
the operator ([3]) and are the same found in Ref. [3(. The 
contribution due to collisions ([5]) reads 

c col = -7 ( ° t Zg ) (8) 
where the auxiliary NxN matrix w is defined by 



Vi+ij + Vi-ij + Vij-i + Vij+i 

«";,: \ v i± ij + v ijTl - 2vij 

Vi-ij-i + v i+ i J+ i - 2vij 



i-j=±l 
i=3 



The matrices u and v are symmetric by construction, 
and one can easily check that z is antisymmetric in the 
stationary state. 



Observables - The most relevant observables are the 
kinetic temperature field Tj = (p?) = va and the local 
energy current Ji , that can be written as the sum of two 
contributions, Ji — jf + Jf, where 



Jf 



oj {qi-iPi 

7 



M) ~ (pU)} 



U Zi-x,i 

7, 



(9) 
(10) 



are the deterministic contribution (due to the springs) 
and the stochastic one (due to collisions), respectively. 

Coordinate change - In the perspective of better un- 
derstanding the differences between the equilibrium and 
non-equilibrium invariant measure, it is convenient to 
choose new variables Qi , Pi in such a way that the equi- 
librium covariance matrix becomes fully diagonal. This 
is accomplished by the linear trasformation 



= aq i+1 - bqi\ P t = p t 



(11) 



with a = [uj 2 + k/2 + {uj 2 k + /c 2 /4) 1 /2]i/2 ; 5 = w 2/ fl _ 
Moreover, in analogy to Eq. ([6]), we introduce a new co- 
variance matrix C whose components U, V, and Z are 
given by U i3 = (QiQj), = {P%Pj), and Zy = {Q l P i ). 

Numerical solution - The stationary solution of Eq. ([7]) 
can be efficiently computed by exploiting the sparsity of 
the corresponding linear problem, as well as the symme- 
tries of the unknowns u, v and z. The elements of C 
can be thereby obtained from transformation (fTTj) . The 
diagonal elements of C correspond to the temperature 
profile. This is illustrated in Fig. QJi, where both Va and 
Uu are plotted versus y = 2i/N — 1 for k = 2 and differ- 
ent system sizes. The linear profile confirms the validity 
of Fourier's law, as expected in the presence of an on- 
site potential 0] ■ All the other C elements are of order 
0(1/N) and proportional to the temperature difference 
AT and are, therefore, scaled accordingly in the other 
panels of Fig. Q] The upper-diagonal of Z is plotted in 
panel (b); because of the antisymmetry of z, Zj_x,i co- 
incides with the deterministic component jf of the flux 
(see Eq. ©). The singular behaviour at the extrema is a 
consequence of Jf vanishing at the boundaries, where the 
particles are not free to move. The behavior of V along 
the upper diagonal is presented in Fig. [TJ:; it basically 
quantifies the 1/N finite-size deviations exhibited by the 
temperature profile. Finally, the dependence of Z along 
the principal antidiagonal [x = (i — j)/N] is shown in 
Fig. [TJi, where one can notice a slow decay of the correla- 
tions. This is analogous to what found in purely stochas- 
tic models [1, HH and is believed to be a generic feature 
of nonequilibrium stationary states [H . 

Principal component analysis - Once the covariance 
matrix is known, it is natural to determine its eigenval- 
ues [y = 1, . . . , 27V) - that we assume to be ordered 
from the largest to the smallest one - and the corre- 
sponding normalized eigenvectors, that are denoted as 

\ . . . , (j>x , V'i^' ■ • • ' ^n^)- Since C is almost diago- 
nal, \^ yields the temperature profile and each of the 
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FIG. 1: (Color online) Some elements of the covariance matrix 
C (properly normalized) for k = 2, ui = 1, 7 = 1, A = 1, 7+ 
I , ' . T- = 0.9 and different chain sizes N = 128,256,512; 
(a) Ui,i = Vi,i = Ti, (b) NZi, i+1 /AT, (c) NVi, i+ i/AT, (d) 
N Zi t N+i-i (antidiagonal) . 



two eigenvector components are localized. A less obvi- 
ous fact is that the eigenvalues come in almost degen- 
erate pairs (the relative difference being much smaller 
than 1/N) and that the <f> and ip components of cor- 
responding eigenvectors are localized around the same 
site m = vjl (see Fig. [2k, where the tp components are 
plotted for v = 100,101). We interpret this by saying 
that any pair (y, v + 1) of eigenvectors identifies a sin- 
gle, localized, degree of freedom at equilibrium with a 
"temperature" A^ . 

The eigenvector spatial structure can be fruitfully il- 
lustrated by plotting pf ] = [4 V) ] 2 + + + 
[V4' y+1 ' ) ] 2 - In Fig- Et>, we see that the vector width docs 
not depend on the system size, which only determines 
the extension of the tails. Moreover, Fig. [2Tj shows also 
that the squared envelope decays as a power-law from its 



localization center, p. 



[£/\i 



hj. This defines 



a typical length £ that can be interpreted as the minimal 
size of the spatial region that is necessary to ensure a 
local thermal equilibrium, a sort of mean free path. This 
interpretation is supported by the fact that £ decreases 
upon increasing the strength 7 of the stochastic process, 
which is the only thermalization mechanism in the bulk. 

The unpinned case - The limit case of vanishing pin- 
ning force, k — 0, is of particular interest. Indeed, 
here the total momentum is conserved and we expect an 
anomalous behaviour, i.e. a diverging finite-size thermal 
conductivity [2, [l(| • As shown in Fig. [3k. the temperature 
profile (solid line) is no more linear. Remarkably, we find 
also that all the nondiagonal elements are in this case of 
order 0(l/\/N) (see the solid lines in Fig. [3Jd,c and d). 
In particular, through Eq. © this implies that the heat 
conductivity diverges as >/N, in agreement with the lin- 
ear response prediction [9( . This observation is supported 
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FIG. 2: (Color online) Localized structure of the eigenvec- 
tors of C. The upper panels correspond to the model and 
parameter values of Fig. 1: (a) V components of the vectors 

v = 100, 101 for N = 200; (b) the squared envelope de- 
fined in the text for N = 400 with 7 = 2 (dashed line) and 
N = 800 with 7 = 0.5, 1, 2 (solid lines from top to bottom), 
(c) P = E^Up'"': FPU model for N = 511, v x = 150, 
V2 = 350 (dotted line), unpinned case (7 = 0.4) for TV = 512, 

vi = 150, v 2 = 350 (solid line) and N = 256, v\ = 75, 
V2 = 175 (dashed line); (d) same data in log-log scale for the 
unpinned case (N — 512). The straight lines have a slope —2 
in (b) and —1.5 in (d). 



by analytical arguments based on multiple-scale expan- 
sion of Eq. ([7]) in the smallness parameter 1/yN 14 1. 
A second remarkable difference with the case k > is 



the exponential decay of the off-diagonal terms, the de- 
cay length scaling as y/N. For instance, in Fig. [3ji (solid 
line) the elements along the antidiagonal of the matrix Z 
are plotted versus x^/N = (i—j) /V~N ; similar curves are 
obtained for the same elements of the matrices U and V. 

On the other hand, the principal component analy- 
sis reveals that the eigenvalues still appear in almost 
degenerate pairs and the eigenvectors are localized (see 
Fig. [2b)) although the tails decay with a slower power-law 
exponent, close to 1.5 (see Fig. [2pl). 

Comparison with a nonlinear model - In order to assess 
the generality of the scenario resulting from the stochas- 
tic model, it is crucial to investigate also a nonlinear de- 
terministic model. We have considered the Fermi-Pasta- 
Ulam (FPU) chain with a purely quartic interparticle po- 
tential (ft+i — (?i) 4 /4 [l5| . We have implemented Eq. |T|), 
by working with stochastic heat baths in the limit of a 
small mass for the reservoir's particles Q. This choice al- 
lows us to use a symplectic integrator and thereby to use 
"large" time steps (10~ 2 ), still obtaining accurate results. 
The covariance matrix is then estimated by time averag- 
ing and is thus affected by the unavoidable statistical 
fluctuations (we have sampled about 5 x 10 7 data points 
separated by one time unit). For a meaningful compari- 
son, the Q variables (fTTj) are chosen to be f2j(g,+i — qi), 
where fli is a renormalized frequency that can be deter- 
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FIG. 3: (Color on line) The same matrix elements as in Fig. 1 
in the absence of an on-site potential k — and for 7 = 
0.4. The three solid lines refer to TV = 128, 256, and 512. 
The dashed lines are obtained by simulation of a chain of 511 
quartic FPU oscillators. 



have studied the fluctuations of the single mode ampli- 
tudes and the correlations between pairs of such modes, 
going beyond the second cumulants. Within the numer- 
ical accuracy, we have not found any deviation from the 
assumption of a purely Gaussian distribution, even in the 
FPU model. We can also conjecture that the Gaussian 
assumption should become exact in the thermodynamic 
limit for model ([I]) . We want to point out that, to our 
knowledge, this is the first case where an explicit nontriv- 
ial, even if approximate, representation of the nonequilib- 
rium invariant measure is given for Hamiltonian models. 
As a further development, the stochastic model should be 
modified to reproduce quantitatively the universal scal- 
ing laws predicted for a generic system [l(| , as well as to 
analyse the corresponding modifications on the structure 
of the measure. 

We acknowledge useful discussions with G. Basile and 
S. Olla. One of us (AP) wishes to thank the Erwin 
Schrodinger Institute for profitable exchanges of ideas. 



mined by imposing that the effective harmonic potential 
energy equals the kinetic energy. The homogeneity of 
the potential suggests setting Q$ = (T^ 4 ; we have found * 
that the best overlap between kinetic and potential en- 
ergy is obtained for £ = 1.4 13(. The results obtained [1] 



for a chain of 511 oscillators are plotted in Fig. [3] (dashed 
lines), where the same y/~N scaling of variables as in the 
unpinned stochastic model is implicitly assumed. Even 
though we should notice that 7 = 0.4 for the stochastic 
model has been selected to yield the best agreement with 
the temperature profile of the FPU model, the similar- 
ity between the deterministic and the stochastic model 
is quite striking and extends to the temperature devia- 
tions as well as to the behavior of the off-diagonal terms. 
The only elements that are significantly different are the 
upper diagonal terms of Z (see Fig. \5^>) . Statistical fluc- 
tuations prevent an accurate analysis of the single eigen- 
vectors. In order to reduce their effects, we have decided 
to average a subset of eigenvectors around their localiza- 
tion center. The results presented in Fig. [2]: (the average, 
denoted by p, is performed over the eigenvectors from 
150 to 350) reveal that also the eigenvectors of the FPU 
model are localized (the side peaks seem to be a finite- 
size effect, as they tend to disappear upon increasing the 
system size) . A quantitative analysis of the decay rate is 
out of question. 

The study of the above models has shown that the key 
features of nonequilibrium steady states are captured by 
principal component analysis and are contained in the 
statement that, the invariant measure can be effectively 
approximated as the product of independent Gaussians 
for the collective, localized mode amplitudes. The vari- 
ances are the local temperatures. 

In order to further explore the validity of this claim, we 
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